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Abstract 

In previous articles we have investigated the firing properties of the 
standard Hodgkin-Huxley (HH) systems of ordinary and partial differen- 
tial equations in response to input currents composed of a drift (mean) 
and additive Gaussian white noise. For certain values of the mean cur- 
rent, as the noise amplitude increased from zero, the firing rate exhibited 
a minimum and this phenomenon was called inverse stochastic resonance 
(ISR). Here we analyse the underlying transitions from a stable equilib- 
rium point to the limit cycle and vice-versa. Focusing on the case of a 
mean input current density /x = 6.8 at which repetitive firing occurs and 
ISR had been found to be pronounced, some of the properties of the corre- 
sponding stable equilibrium point are found. A linearized approximation 
around this point has oscillatory solutions from whose maxima spikes tend 
to occur. A one dimensional diffusion is also constructed for small noise 
based on the correlations between the pairs of HH variables and the small 
magnitudes of the fluctuations in two of them. Properties of the basin 
of attraction of the limit cycle (spike) are investigated heuristically and 
also the nature of distribution of spikes at very small noise correspond- 
ing to trajectories which do not ever enter the basin of attraction of the 
equilibrium point. Long term trials of duration 500000 ms are carried out 
for values of the noise parameter a from to 2.0, with results appearing 
in Section 3. The graph of mean spike count versus a is divided into 4 
regions R\,...,Ra, where R3 contains the minimum associated with ISR. 
In Ri noise has practically no effect until a critical value of a — a ci is 
reached. At a larger critical value a — a C2 , the probability of transitions 
from the basin of attraction of the equilibrium point to that of the limit 
cycle becomes greater than zero and the spike rate thereafter increases 
with increasing a. The quantitative scheme underlying the ISR curve 
is outlined in terms of exit time random variables and illustrated dia- 
grammatically. In the final subsection 3.4, several statistical properties of 
the main random variables associated with long term spiking activity are 
given, including distributions of exit times from the two relevant basins 
of attraction and the interspike interval. 

Short Title: Hodgkin-Huxley 

Keywords and Phrases: Hodgkin-Huxley equations, noise 

1 Introduction 

The Hodgkin-Huxley [1] systems of ordinary and partial differential equations, 
based on the electrophysiology of the squid giant axon, are the cornerstone of 
mathematical models of single neurons as well as several types of cardiac cells. 
Recent such studies include those of Komendantov et al. [2] for hypothathalamic 
magnocellular neuroendocrine cells, Saarinen et al. [3] for cerebellar granule 
cells, Williams et al. [4] for ventricular myocytes, Kameneva et al. [5] for retinal 
ganglion cells and Drion et al. [6] for dopaminergic neurons. Many of these 
computational cell models contain 10 or more components as the important 
roles of many different ion channels have been discovered since the appearance 
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of the HH model. Analysis of such higher-dimensional models is very complex 
as there may be 50 or more parameters in distinction to the relatively few in 
the 4-component Hodgkin-Huxley system. Since the latter does in fact embrace 
some of the basic firing properties of neurons in general, there has naturally been 
a large number of analyses and computational studies of the HH systems. These 
include both deterministic (for example, [7, 8, 9, 10, 11, 12]) and stochastic (for 
example, [13, 14, 15, 16, 17]) modeling. 

In recent articles [19, 20, 21] we have explored the effects of both additive 
Gaussian white noise and conductance noise on repetitive firing in the Hodgkin- 
Huxley system. In the additive noise case for the ordinary differential equation 
(ODE) model, when the mean input current density /i is not far above the 
threshold of 6.4 fiA/cm 2 for repetitive firing, the number of spikes in the first 
500 or 1000 ms was found to undergo a pronounced minimum (ISR) as the noise 
level a increased from zero [19]. The minimum occurred around a = 0.35. Guo 
[22] has recently found similar results for the HH ODE system with colored 
(Ornstein-Uhlenbeck process) noise. Similar results were found for the partial 
differential equation (PDE) system [20, 21] where the spatial distribution of 
the noise was also an important factor, which led to the disinction between the 
effects of noise on the instigation and propagation of spikes. 



1.1 Model description 

In this article we restrict attention to the HH ODE system, which corresponds 
to a uniformly polarized or "space-clamped" neuron. The system of stochastic 
differential equations was given in our previous articles but are repeated here 
for completeness and notation: 

dV = i \fi + g K n\V K -V) + g Na m 3 h(V Na -V)+ g L (V L - V)]dt + adW (1) 

and for the auxiliary variables 

dn = [a n (l — n) — j3 n n]dt (2) 

dm = [a m {l — m) — /3 m m]dt (3) 

dh = [a h (l - h) - /3 h h]dt (4) 

where C is the membrane capacitance per unit area, /i, which may depend on t, 
is the mean input current density, g Kl g Na and gL are the maximal (constant) 
potassium, sodium and leak conductances per unit area with corresponding 
equilibrium potentials Vk, Vjva> an d Vi, respectively. The noise enters as the 
derivative of a standard Wiener process W and has amplitude a. The auxiliary 
variables are n(t), the potassium activation, m(t), the sodium activation and 
h(t), the sodium inactivation. The coefficients in the differential equations for 
the auxiliary variables as functions of depolarization are 

aniY) = 100[e( 10 -^)/ 10 - 1] (5) 
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Pn(V) = \e- v '™ (6) 

n[ ' ~ 10[e( 25 - y )/ 10 - 1] { ) 

P m {V) = Ae~ v ' w (8) 



2 Stable equilibrium point and limit cycle 

When fi is above the critical value /i Cl for repetitive firing (saddle-node bifur- 
cation) and smaller than the value ^ C2 at which there is a subcritical Hopf 
bifurcation, there are two attractors consisting of a limit-cycle (action potential 
trajectory) and a stable equilibrium point. 

2.1 Stable equilibrium point x* 

Let us denote the random vector (V, n,m, h) by X = (Xi, X 2 , X 3} X4) and 
rewrite the system of equations as 

dXi = F!(X)dt + adW (11) 

dX k = F k (X)dt (12) 

where k = 2, 3, 4. Because a mean input current density of fj, =6.8 /zA/cm 2 was 
found to exhibit a pronounced minimum in firing as a increased, we will focus 
on this value of /x unless stated otherwise. 
Setting, with zero noise, 

F(x*) = (13) 

yields equilibrium points for the deterministic Hodgkin-Huxley system. 
The Jacobian matrix at the equilibrium point is defined as 

j(x * )= {SL' u=1 '-' 4 - (i4) 
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Figure 1: Eigenvalues in the complex A-plane for three values of \x. The values 
of the real eigenvalue near A = —0.13 are indistinguishable in this diagram. 

Figure 1 shows the eigenvalues in the complex lambda plane for values of 
ji = 5, below the critical value /i Cl , /i = 7.5, between the critical values /j, Ci , 
and /x C2 , and fi = 10 which is above /i C2 . The nature of the critical point at 
the various values of fi can be readily seen. For the chosen value of /i = 6.8, 
numerical evaluation gives an equilibrium point at 



at which the components of F are 

(-0.00000003084, 0.000024246, -0.0000013509, 0.0000048941) (16) 

Numerical evaluation of the eigenvalues of J(x*) gives Ai = —4.641, the complex 
conjugate pair A2,3 = —0.630 ± 0.548i and A 4 = —0.1323. Hence x* is an 
asymptotically stable spiral point. 

2.2 The linear approximation 

The system of stochastic ordinary differential equations for the process X ob- 
tained by linearizing about x* is 



x* = (4.0536,0.38107,0.084327,0.45129) 



(15) 




(17) 
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of, 



dXj =J2g^Xkdt,j = 2,3,4, 



fe=i 



(18) 



where the partial derivatives are evaluated at x* . The Jacobian is found nu- 
merically to be 



10 2 * 



-0.010891 -1.2764 1.2710 0.079467 

0.000030551 -0.0019202 

0.00032794 -0.034876 

-0.000044773 -0.0012664 



(19) 



Note that the system (17)-(18) is linear and so X is a Gaussian process. 
Using the theory in Rodriguez and Tuckwell [23], in the absence of an imposed 
spiking threshold, the exact distribution of X can be found at any t. 

Between spikes the linearized system (17)-(18) is expected to provide a rea- 
sonable approximation to the fully nonlinear system (l)-(4). This is clearly 
demonstrated by the two sets of sample paths shown in Figure 2. With input 
parameters /i = 6.8, a — 0.6, a time segment of length about 60 ms gave the 
sample path for V (Xi) shown in the top part of the Figure. With the same 
path for the Wiener process (or the white noise), the sample path for the voltage 
X\ in the linearized system, with the same initial values, is seen in the lower 
part of Figure 2, to mimic closely that in the upper part. However, there is one 
very striking difference betwen the two paths as at about t = 59 ms the voltage 
in the linear system attains a local maximum and then decreases quite rapidly, 
continuing in an oscillatory fashion. In contrast, at about the same t and volt- 
age values, a spike arises in the nonlinear system. However, these paths indicate 
that the time of spiking in the nonlinear system can be well approximated as 
the time at which the voltage in the linear sytem first attains a threshold value. 
Such a first passage time can be determined from the usual first exit time theory 
for diffusion processes. 
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Figure 2: In the top part, voltage is plotted versus time for a sample path in 
the nonlinear full Hodgkin. Huxley system with the parameters shown. A spike 
forms near the end of the record. In the lower record is shown the voltage path 
for the system linearized about the stable equilibrium point. The Wiener path 
is the same in both records to enable a comparison to be made. 



2.3 Numerical examples 

Without an imposed threshold condition, no spikes are possible in the linearized 
stochastic system described by (17) and (18) because there is only a stable equi- 
librium point about which trajectories fluctuate. This is of course in distinction 
to the Hodgkin-Huxley system (with suitable input parameters) where trajec- 
tories may, if the fluctuations are large enough, give rise to spikes around the 
limit cycle. It is of interest to examine some statistical properties of the original 
process X in nonspiking periods. An example of paths with fi = 6.8 and a = 0.1 
over a 50 ms nonspiking period is shown in Figure 2. The basic statistical prop- 
erties of the components over a 500 ms time period are given in Table 1. In 
Table 2 are given the corresponding correlation coefficients. 
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Table 1: Statistics of variables during non-spiking 
period, fi — 6.8, a = 0.1 





Mean 


Max 


Min 


St Dev 


Coef Var 


V 


4.049 


4.446 


3.658 


0.1337 


0.0330 


n 


0.3811 


0.3828 


0.3795 


6.61e-4 


0.0017 


m 


0.0843 


0.0874 


0.0811 


0.0012 


0.0143 


h 


0.4515 


0.4542 


0.4488 


0.001 


0.0023 



Table 2: Correlation coefficients during non-spiking 
period, fi — 6.8, a = 0.1 





V 


n 


m 


h 


V 


1.0000 


0.3068 


0.9562 


-0.2137 


n 


0.3068 


1.0000 


0.4649 


-0.9894 


m 


0.9562 


0.4649 


1.0000 


-0.3699 


h 


-0.2137 


-0.9894 


-0.3699 


1.0000 



0.382 

n 

0.38 



0.086 
m 0084 
0.082 



0.45 
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Figure 3: Sample paths for the 4 components of X over a 50 ms period during 
which there were no spikes. Input parameters fi = 6.8 and a = 0.1. 

Tables 3 and 4 give the statistical properties during a 500 ms nonspiking 
time period when the noise level is a = 0.6. 
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Table 3: Statistics of variables during non-spiking 
period, fi — 6.8, a = 0.6 





Mean 


Max 


Min 


St Dev 


Coef Var 


V 


4.149 


8.962 


-0.206 


1.301 


0.314 


n 


0.383 


0.419 


0.368 


0.0077 


0.020 


m 


0.0859 


0.1393 


0.0538 


0.0126 


0.147 


h 


0.4471 


0.4668 


0.3883 


0.0122 


0.0274 



Table 4: Correlation coefficients, [i = 6.8, a = 0.6 





V 


n 


m 


h 


V 


1.0000 


0.3505 


0.9690 


-0.2408 


n 


0.3505 


1.0000 


0.5284 


-0.9854 


m 


0.9690 


0.5284 


1.0000 


-0.4272 


h 


-0.2408 


-0.9854 


-0.4272 


1.0000 



There are two very noticeable features of the sample paths in the nonspiking 
periods examined. Firstly, the oscillatory character of the paths, which is trace- 
able to the eigenvalues of the Jacobian at the equilibrium point. Secondly, the 
strong positive correlation between V and m and the strong negative correlation 
between n and h. 




Figure 4: Showing how single, double and multiple spikes arise from oscillations 
in the nonlinear full Hodgkin-Huxley system with noise, /i = 6.8 
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The oscillatory nature of the paths can never be captured in the standard 
integrate and fire models nor the leaky integrate and fire models [24]. When 
the noise is sufficiently large to make for fairly frequent spiking, the times of 
spiking must tend to arise at the maxima in the oscillations of V. This is seen 
dramatically in Figure 3 and in Figure 2 of the previous subsection. 

2.4 A one-dimensional diffusion 

Examination of the statistical properties of the process as given in Tables 1 
and 2 shows that for small noise the coefficients of variation of n and h are an 
order of magnitude smaller than those for V and m and that the correlation 
coefficient of V and m is close to unity. These observations make it reasonable 
to consider a 1-dimensional approximation V to the HH system in which a 
constant n w n and a constant h w h, with m = kV, where k is another 
constant. These approximations lead, on putting C = 1 in (1), to the following 
stochastic differential equation for V, 

dV = {n + Cl - c 2 V + c 3 V 3 - c 4 V A )dt + adW, (20) 

where 

ci = g K n 4 V K + g L V L (21) 

C2 = Qk 1 ^ + 9l (22) 

c 3 = g Na k 3 hV Na (23) 

c 4 - g Na k 3 h. (24) 

Standard theory gives for such a diffusion that the mean exit time from a value 
x € (a, b) to outside this interval satisfies the ordinary differential equation 

(a 2 /2)M" + (/i + ci - c x + c 3 x 3 - 4)M' = -l,x€ (a, 6), (25) 

where primes denote differentiation, along with suitable boundary conditions. 
Preliminary investigations of the validity of this approximation were made using 
values for /j, = 6.8 and a = 0.1 (sec Tables 1 and 2) but a detailed study will be 
reported in a subsequent article. 

2.5 The limit cycle and its basin of attraction 

With constant input current density in the interval \p Cl , jj, C2 ) , repetitive spiking 
may occur with a fixed period. For fj, = 6.8 and no noise the period is about 
17.65 ms. The limit cycle is in 4-space but we here show in the upper part of 
Figure 2 the projection of the limit cycle obtained by plotting n versus V. In 
the lower part of the figure is shown the position of the stable equilibrium point, 
here designated R, for the same value of /i. It can be seen that the limit cycle 
approaches quite close to R. In a previous article in which the moment method 
was used to explore the effects of noise on HH spiking (Tuckwell and Jost, 2009), 
we heuristically estimated the basin of attraction of the stable equilibrium point. 
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Part of the basin of attraction of the limit cycle can be numerically estimated 
by taking the union of all stochastic paths which do not collapse to the stable 
equilibrium point. As explained later, there are ranges of values of a where the 
stochastic paths have this property and Figure 6 shows a sample of such paths 
for various such a. Here paths were taken over 1000 msec involving about 50 
spikes for values of a ranging from 0.05 to 0.25. This picture gives an idea of 
how far off the deterministic limit cycle a path may wander without entering 
the basin of attraction of the stable equilibrium point. The variability of paths 
is relatively small. 




0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

V(mV) 




n 



Figure 5: In the top part, voltage is plotted versus potassium activation vari- 
ablefor the deterministic path of repetitive spikes, depicting the limit cycle for 
the HH ODE system with /i = 6.8 > fi Cl . In the lower part the limit cycle is 
magnified in the vicinity of the stable rest point. 
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Figure 6: The union of many stochastic paths for values of the noise parameter 
cr for which paths did not collapse into the stable point in a time interval of at 
least 1000 ms. 

Another way to see the limited variability of the times to complete a spike 
orbit (limit cycle) is displayed in Figure 7. Here spikes were observed during 
repetitive firing during periods in which no transitions to the basin of attraction 
of the rest point occurred, with noise levels from a — to a — 0.4. The statistical 
properties of the interspike interval (ISI) are shown in the Figure. The most 
salient feature is that the mean ISI is practically constant (blue triangles) as the 
noise varies, staying in the interval [17.58, 17.82] for these values of a. Naturally 
the standard deviation of the ISI increases (green squares) as a increases, the 
maximum ISI increases /red circles) and the minimum ISI decreases (black 
diamonds). 
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Figure 7: Some statistical measures of the ISI during repetitive spiking at various 
noise levels with /i = 6.8. For a = 0, there is no variability. The mean ISI is 
shown with (blue) triangles, the (green) squares denote the mean + 3 standard 
deviations, the (red) circles denote the maxima and the (black) diamonds denote 
the minima. Until a is about 0.07, there are no cessations of spiking up to 500000 
ms. 

Examples of distributions of the ISI during repetitive spiking at small noise 
levels are shown in Figure 8. In the top panel is an ISI histogram for a noise level 
(<r = 0.07) at which there is apparently no cessation of spiking over extremely 
long (infinite?) time periods with 28429 ISIs in 500000 ms. The distribution is 
roughly Gaussian and has a mean of 17.59 ms and a standard deviation of 0.221 
ms. For the second histogram shown, a = 0.085, which is just greater than the 
critical value at which spiking may stop in a finite time. The number of spikes 
is only 4397 and the mean and standard deviation of the ISI are 17.60 ms and 
0.276 ms, respectively. The distribution of the ISI is practically the same in 
both cases. 
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Figure 8: Histograms of the ISI during repetitive spiking at two small noise 
levels. Top, a = 0.070; bottom, a — 0.085. For a = 0.07 there are no cessations 
of spiking up to 500000 ms whereas for a — 0.085 firing stops after about 4400 
spikes. 



3 Results on spiking for the HH system and in- 
verse stochastic resonance 

In the following it is assumed that [i is such that repetitive spiking does in 
fact occur. In relation to the stochastic paths for the HH system of stochastic 
equations we define the following two random variables. 

Firstly, the exit time of the process to escape from the basin of attraction 
Bl of the limit cycle L to that Br of the rest point R is by 

T L ^ R {x ),x Q eB L (26) 

where (xq) = (Vo, no, mo, ho) is an initial point. Secondly the exit time for the 
process to escape from the basin of attraction Br of R to that of the limit cycle 
is 

T R ^ L (x ),x Q eB R . (27) 
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Using the standard theory for diffusion processes, Kolmogorov second order par- 
tial differential equations for the moments and distributions of these quantities 
as a function of initial values were described in our previous work (Tuckwell et 
al., 2009). However, any attempt to solve these equations analytically or even 
by numerical methods for PDEs, apart from being a formidable task, requires 
an exact or even approximate knowledge of the two basins of attraction which is 
unfortunately not presently available. In fact the probabilistic nature of Bl and 
Br is not completely understood because it seems that for small enough noise, 
at least for some values of /x, escape from Bl may be impossible, not being 
observed in extremely long time periods and similarly, for particular ranges of 
a, for escape from Br. 

3.1 Long term trials and data collected 

We are interested in obtaining samples of meaningful sizes for the above random 
variables and determining their properties, mainly as a function of a . Since some 
of the exits from one basin of attraction to the other are, for some values of <r, 
very rare events the simulations of solutions of the full stochastic HH system 
have to be performed over very long time periods. These were chosen to be 
500000 ms or about 8 minutes and 20 seconds with a timestep of 0.065 ms. 
Results for the 500000 ms interval were obtained as the union of 100 trials of 
length 5000 ms, where the final values of V,n,m,h from any trial were used 
as the initial values for the next trial. Then the data on t, V, n, m, h for the 
100 trials were concatenated to give 5 single vectors each with over 7.5 million 
elements. The records for V were analyzed to determine the times of occurrence 
of spikes and from these the interspike intervals were determined. In total, 50 
trials of length 500000 were simulated. 

3.2 All spikes 

The numbers of spikes were recorded in each of 50 trials of length 500000 ms 
for 40 values of a. The mean number of spikes, denoted by E[NSP], is plotted 
against a in Figure 9. In the top panel all results are included, whereas in the 
lower left panel detail is shown for < a < 0.1 and in the lower right panel, for 
0.1 < a < 0.5. 

The general form of the plot of E[NSP] versus a is similar to that in Figure 
5 of Tuckwell et al. (2009) where the total time period was much less at 1000 
ms. Thus, Figure 9 exhibits the phenomenon of inverse stochastic resonance, 
which refers to a firing rate which, as noise level increases, at first declines to 
a minimum and then becomes greater. However, Figure 9 shows more detail 
and reveals 4 distinct regimes marked R\,...,Ri. In region R\, which extends 
from the deterministic setting of a = to very close to a = 0.07, the noise 
is almost without effect and the number of spikes is always close to that for 
the zero noise, 28431. Region i? 2 is characterized by a rapidly falling number 
of spikes as a increases from 0.07 to 0.14. At the latter value of a the mean 
number of spikes is 104.8 which is about 0.4% of the maximum value. In region 
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i?3) from a = 0.15 to a = 0.35, E[NSP] stays below 100 with a broad minimum 
of about 9.5 spikes (0.03% of the maximum) around a = 0.295 to 0.300. By 
a = 0.375, the beginning of region R4, the mean number of spikes has reached 
over 120 and then increases sharply and eventually at a slower rate to reach 
25883 around a = 2. 



E[NSP] 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 



E[NSP] 




0.05 0.1 0.1 0.2 0.3 0.4 0.5 



Figure 9: The dependence on the noise parameter a of the mean of the total 
number of spikes E[NSP] in a time interval of length 500000 ms in the nonlinear 
HH system. Here /i = 6.8 with 50 trials at each point. The values of a arc 
divided into 4 regimes designated R\ to R4. In the top panel, results are given 
for all values of a. In the lower left panel the detail of R\ and the start of R2 
are shown. In the lower right panel, the detail of R3 and the start of R4 are 
shown. 



3.3 Underlying scheme for ISR 

The following observations constitute a basis for ISR which occurs at some values 
of fi in the HH system described by Equations (l)-(4). We define two critical 
values a Cl and a C2 of the noise parameter er, with < a Cl < er C2 < 00. er Cl and 
<7 C2 depend on /1, and are only relevant above the critical value for repetitive 
firing and possibly for selected initial values of the process (see Tuckwell et al., 
2009). With reference to the two random variables defined by (26) and (27), 
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but without any specific values of Xq, we have the following, as schematized 
diagrammatically in Figure 10. 
< er < er Cl . 

Repetitive spiking continues, presumably indefinitely. 

The probability of a transition from Bl to Br at any time is zero. 

The expectation of Tl->r(xo), xq e Bl, is infinite. 

cr a < f < f C2 . 

Repetitive spiking ceases at some time ii < oo. 

The probability of a transition from Bl to Br is greater than zero and increases 
with noise level. 

The expectation of Tl->r(xo), x € Bl, is finite and and tends to decrease with 
increasing a. 

After a transition from Bl to Br, the probability of the reverse transition from 
Br to Bl is zero. 

The expectation of Tr^l( x o), x o € Br, is infinite. 

After the cessation of repetitive spiking, there is no more spiking. 

cr > rr, , . 

After a transition from Bl to -B/j, the probability of a transition from to 

Bl is greater than zero and tends to increase as a incresases. 

The expectation of T r _+ l (xq), Xq e Br, is finite and tends to decrease as cr 

increases. 

After a transition back from Br to Bl, the reverse transition occurs with non- 
zero pobability, followed again by a transition from Bl to Br and so forth, 
indefinitely 

Epochs of spiking and nonspiking alternate with variable durations which tend 
to become shorter as a increases. 

Figure 10 illustrates some of these aspects. The bottom panel shows part of 
the actual E[NSP] curve (normalized) from Figure 9 and the approximate values 
of a Cl and <r C2 are obtained from this curve. In the top panel, the probabilities 
of transitions from Bl to Br and Br to Bl are indicated as commencing to 
increase from zero at the two critical values of a and rising to unity, although 
this final value depends on the initial value x which is not fixed throughout. 
In the middle panel, descriptions of the expectations of the exit times from Bl 
and Br are given. 
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Figure 10: A schematic representation for the transitions underlying ISR. There 
are shown two critical values of the noise parameter a and the probabilities of 
transitions between the basins of attraction and Br of the limit cycle and 
rest point are sketched in the top panel for the various ranges of values of a. In 
the bottom panel is shown the normalized expected number of spikes obtained 
by simulation over 500000 ms for /i = 6.8, a value of fi at which ISR had been 
found to be pronounced. 



The scheme given in Figure 10 is in accordance with the results for sample 
paths of V for simulated spike trains shown in Figure 11. Here the time periods 
are all 500000 ms but in the top two records time is only shown to 2000 ms. The 
very small noise case gives apparently indefinite repetitive spiking. Somewhat 
larger noise leads to an initial isolated burst followed by silence, possibly for 
an infinite time. Increasing the noise further leads to occasional bursting and 
eventually the frequent bursting case occurs. 
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Figure 11: Illustrating the 4 basic patterns of spiking activity. Top left: a — 
0.05. Incessant spiking, at least to 500000 ms. Top right: a — 0.25. One initial 
burst with no further spikes, at least to 500000 ms. Bottom left: a — 0.34. 
Occasional bursts with separations of order 100000 ms. Bottom right: a = 0.45. 
Frequent bursts. 



3.4 Some statistics of T l ^r and Tr^ l 

In this subsection results on certain statistical aspects of the underlying random 
variables T^r and Tr-^Zi defined by (26) and (27), obtained from the long- 
term simulation of the HH system (l)-(4) are given. 



3.4.1 Escape time from Bl 

During repetitive spiking the trajectories of the process do not deviate much 
from the deterministic limit cycle as illustrated in Figure 6. It is assumed that 
there is a well-defined but unkown set in (V, n, m, /i)-space, denoted by B^, 
containing the deterministic repetitive spiking trajectory, which is called the 
basin of attraction of the limit cycle. This implies that if the process started in 
Bl, then, with no noise, the path would approach the limit cycle. With noise, 
trajectories may escape from Bl whereupon they enter the basin of attraction 
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Br of the stable equilibrium point. The nature of Bl is unknown and it is 
not clear that it is a regular open set because if it was it is likely that escape 
from it would occur eventually with probability one no matter how small the 
noise level. Thus, even though it seems that for a < cr Cl , P(L —> R) = and 
E[Ti,_># = oo], it may be the case that P(L — > R) is so small for small enough a 
that the event of this escape is unlikely to be observed in the course of feasible 
simulations. It is a remaining mathematical challenge to ascertain the veracity 
of these remarks. 

Notwithstanding these uncertainties, the mean value of E[Ti_j.^] of the exit 
time from B^ was estimated from sample paths and the results are shown for 
a > 0.2 in Figure 12. For these noise levels the mean exit time declines rapidly 
from very large values to a minimum of about 57 ms at a — 1.25. Thereafter, 
E[Tl^.r] seems to increase slightly to about 72 ms at a — 2. 
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Figure 12: The dependence on the noise parameter a of the expectation E[7l^/?] 
of the random variable which is the time of exit of the process from the basin 
of attraction of the limit cycle to that of the stable rest point. 50 trials at each 
point with /i = 6.8. For values of a less than a critical value <r Cl « 0.07 the 
value of this expectation is apparently infinite. Extremely large values which 
occurred for a Cl < a < 0.2 are not shown. 
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3.4.2 Escape time from Br 

The basin of attraction of the stable rest point is also unknown exactly. Figures 
13 and 14 show estimates of the mean and standard deviation of Tr_j,£ over 
various ranges of values of a. Again, it is not known with certainty, but it 
appears that the probability of escape from Br to Bl is zero (or extremely 
close to zero) until a > a C2 . In Figure 13, the mean is shown dropping from 
values of order 300000 ms at a = 0.25 to eventually reach values about 30 ms 
at <T = 2. Figure 14 shows corresponding results for the standard deviation, but 
only for a > 0.35 because sample sizes were too small for lesser values of the 
noise parameter. In the lower panel of Figure 14 is shown the dependence of the 
coefficient of variation of Tr^,l on noise level. It can be seen that there is an 
initial increase in this quantity until a maximum is attained at about a = 0.5, 
whereupon it declines monotonically. 




c 



Figure 13: The dependence on the noise parameter a of the expectation E[T#_^£] 
of the random variable which is the time of exit of the process from the basin 
of attraction of the stable rest point to that of the limit cycle. Here fi = 6.8. 50 
trials at each point. For values of a less than some value just less than 0.25 the 
value of this expectation is apparently infinite. 
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Figure 14: Top two panels. The dependence on the noise parameter a of the 
standard deviation o-T R ^ L ,oi the random variable which is the time of exit of 
the process from the basin of attraction of the stable rest point to that of the 
limit cycle. Here fj, = 6.8. 50 trials at each point. Bottom panel. The coefficient 
of variation of the random variable Tr^l as a function of a. 



3.4.3 Distributions 

There are three main random variables of interest in relation to the empirical 
spike trains obtained in this study. These are the previously defined exit times 
Tl^r and Tr^l and in addition the general ISI for the whole train. Samples for 
the latter random variable can be viewed as the union of those of the previous 
two. Examples of histograms of these random variables are shown in Figure 15. 
In all cases the results for 50 trials of duration 500000 ms are combined. In 
Figure 15A is given a histogram for all ISIs for a = 0.2. This occurs in region 
i?3 of Figure 9, where there are very few spikes, all with short ISIs as there are 
no returns to Bl fron Br. The histogram is of the same nature as those in 
Figure 8. See also Figure 7. 

In Figure 15C the raw histogram is shown for all ISIs with a = 0.350. There 
is a preponderance of small intervals and then an exponential-type distribution 
of longer (> 21.5 ms) intervals as depicted in Figure 15E. The tail of the expo- 



22 



nential distribution is very long and extends out to about 500000 ms which is 
the limit in these simulations. The distribution depicted in Figure 15E is close 
to the actual distribution of Tr^l for this value of a. Similar remarks apply to 
the pairs of Figures 15B and 15D (a = 0.375) and 15F and 15G (a = 0.5). In 
the final plot of Figure 15H, the number of spikes has become enormous and the 
majority of ISIs are less than 50 ms as expected for a — 1 which is in region R4. 
of Figure 9. In the histograms of Figure 15 C, B, F and H, there are small bin 
counts out to very large times and these are not visible compared to the large 
bin count at short intervals. However, they are visible at relatively small values 
in the truncated histograms. 
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Figure 15: Examples of distributions (histograms) of ISIs and exit times from 
Br based on 50 trials (data pooled) of length 500000 ms. Red histograms, raw 
data. Blue histograms, ISIs truncated at 21.5 ms. A. Raw data for a = 0.2. 
B. Raw data for a = 0.375. C. Raw data for a — 0.350. D. Truncated data 
for a = 0.375. E. Truncated data for a = 0.350. F. Raw data for a = 0.5. G. 
Truncated data for a — 0.5. H. Raw data for a = 1.0. 

To complete the picture for the distributions of the three key random vari- 
ables, Figure 16 shows histograms of mean numbers of spikes per burst for two 
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values, 0.4 and 1.0, of a. If these arc multiplied by the mean ISI within bursts, 
which is about 17.6 iris, the mean times spent in the basin of attraction of 
the limit cycle {Bl) and hence of Tl^r are estimated. For smaller a some 
very large means were obtained (not shown). Noticeable in Figure 16 is the 
smaller magnitude and small variability when a = 1.0 compared to a = 0.4. 
The distributions in both cases shown are approximately Gaussian which may 
be compared with the exponential-types shown in Figures 15 D, E and G. 
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Figure 16: Histograms of mean numbers of spikes per burst, which indicate 
magnitudes of exit times from Bl to Br. based on 50 trials of length 500000. 



4 Discussion 

The inhibitory effects of noise on repetitive spiking in squid axon, on which the 
HH system is based, were well documented experimentally by Paydarfar et al. 
[25] . The first theoretical evidence of ISR in the full HH sytem was obtained for 
fi = 5 in the non-repetitive spiking mode [26] (Figure 2), as it was found that 
there was a maximum in the mean ISI at small values of a. Subsequently ISR has 
been demonstrated for repetitive spiking in both the HH system of ODEs and 
PDEs [18, 20, 21, 22]. In addition to using the standard initial conditions and 
additive noise, in [18] the phenomenon was confirmed with respect to random 



24 



initial conditions and conductance-based input. Some theory of ISI was also 
outlined [18] in terms of bifurcations in the HH sytem [7, 10, 27] and of the exit 
times [28] from the basins of attraction of a stable equilibrium point and a limit 
cycle. However the simulations were over relatively short time periods and the 
statistical details of the exit times were not explored in detail. 

In this article we have first explored the properties of the stable equilibrium 
point in depth. Linearizing the stochastic HH system about that point gave an 
approximate sytem of stochastic differential equations whose oscillatory solu- 
tions were found to mimic those of the full system in nonspiking periods. The 
oscillations are evidently an integral part of firing in the HH system with the 
considered parameters, because spikes did tend to emerge at the maxima. The 
distribution of ISIs with trajectories not departing greatly from the limit cycle 
for very small noise was estimated and is, as perhaps expected, Gaussian-like. 

Long-term simulations, to 500000 ms were performed with /i = 6.8 for many 
values of the noise parameter a. Four ranges of values of a were distinguished, 
based on mean total numbers of spikes. These regions were denoted R\, i?4 
in Figure 9. Two critical values of a were also apparent, denoted by < <r Cl < 
<r C2 < oo. A detailed discussion of the mechanisms underlying ISR in terms of 
these critical values was given in Section 3.3. It is not certain whether escape 
from the basin of attraction of the limit cycle can ever occur when a < <r Cl or 
whether escape from the basin of attraction from the basin of attraction of the 
stable equilibrium point can ever occur when a Cl < a < cr C2 , but these events 
were never observed in 50 trials of length 500000 ms with several values of a. 
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